LoadXsense

Data visualisation and management

This is the code used to visualise and manage the output of the Xsense dot IMU data. The source code can be found on my github

Load the data

Lets start by defining a function to correctly load measurements:

Note that IMU’s do not start and stop measureing at the exact same time; even after synchronization the amount of elements per IMU (the length of measurement) differs. In my implementation of temporal relaignment I assumed that the time in SampletimeFine was synchronized, and excluded first or last elements accordingly to ensure dataframes are of equal length. This eases calculation since R prefers to calculate over lists of equal length.

# Clean workspace and load dependencies
rm(list = ls())
library(tidyr)
library(dplyr)
library(plotly)
library(ggplot2)
library(knitr)
library(lubridate)

# Here i defined what test subject refers to what directorty
#! should probably think about PID here later on. 
# format: name, hz, skipheading

pp_info <- data.frame (
  #       file, hz, skiprow
  bart = c("../../Logs/old/20240429_163145_bart/", 60, 7),
  other = c("../../Logs/new/20240502_192335/", 60, 10)
)

LoadXsenseData <- function(nameofpp) {
  dir <- nameofpp[1]
  hz <- as.numeric(nameofpp[2])
  skiprow <- as.numeric(nameofpp[3])
  #dir <- bart
  #hz <- 60

  files <- list.files(path = dir, full.names = TRUE)
  data <- list()

  # Read CSV files of each directory
  for (i in seq_along(files)) {
    data[[i]] <- read.csv(files[i], header = TRUE, skip = skiprow)
  }

  # Ensure all dataframes have the same number of rows
  min_rows <- min(sapply(data, nrow))
  data <- lapply(data, function(df) {
    df <- df[1:min_rows, , drop = FALSE]
    return(df)
  })

  # Adjust time
  for (i in seq_along(data)) {
    rows <- nrow(data[[i]])
    data[[i]]$TimeS <- ((1/hz) * (1:rows))
  }

  # Initialize toreturn data frame with time column
  toreturn <- data.frame(time = data[[1]]$TimeS)

  # Calculate absolute values
  for (i in 1:length(data)) {
    if ("FreeAcc_X" %in% names(data[[i]])) {
      col_name <- paste0("FreeAcc_abs", i)
      toreturn[[col_name]] <- sqrt(data[[i]]$FreeAcc_X^2 + data[[i]]$FreeAcc_Y^2 + data[[i]]$FreeAcc_Z^2)
    }
    if ("Acc_X" %in% names(data[[i]])) {
      col_name <- paste0("A_abs", i)
      toreturn[[col_name]] <- sqrt(data[[i]]$Acc_X^2 + data[[i]]$Acc_Y^2 + data[[i]]$Acc_Z^2)
    }
    if ("Gyr_X" %in% names(data[[i]])) {
      col_name <- paste0("Gyr_abs", i)
      toreturn[[col_name]] <- sqrt(data[[i]]$Gyr_X^2 + data[[i]]$Gyr_Y^2 + data[[i]]$Gyr_Z^2)
    }
  }

  # Order the attributes of the dataframe
  toreturn_sorted <- toreturn[, sort(names(toreturn))]

  return(toreturn_sorted)
}

Visualization

lets also define some functions to visualize the data

# Some functions to visualize the acceleration and the Gyr

plot_a <- function(df) {
  plot <- plot_ly(df, x = ~time, y = ~A_abs1, name = "marker1", type = "scatter", mode = "lines") %>%
          add_trace(y = ~A_abs2, name = "marker 2") %>%
          add_trace(y = ~A_abs3, name = "marker 3") %>%
          add_trace(y = ~A_abs4, name = "marker 4") %>%
          add_trace(y = ~A_abs5, name = "marker 5") %>%
          layout(title = "Absolute accelerations",
                 xaxis = list(title = "Time"),
                 yaxis = list(title = "A_abs Values"))
  return(plot)
}

plot_gyr <- function(df) {
  plot <- plot_ly(df, x = ~time, y = ~Gyr_abs1, name = "marker1", type = "scatter", mode = "lines") %>%
          add_trace(y = ~Gyr_abs2, name = "marker 2") %>%
          add_trace(y = ~Gyr_abs3, name = "marker 3") %>%
          add_trace(y = ~Gyr_abs4, name = "marker 4") %>%
          add_trace(y = ~Gyr_abs5, name = "marker 5") %>%
          layout(title = "Absolute Gyr",
                 xaxis = list(title = "Time"),
                 yaxis = list(title = "Gyr Values"))
  return(plot)
}

#! Maybe include a plotting function that takes the dataframe and the attribute to plot, assuming 5 markers?

Initial test

Lets load some data and see how it looks:

note: to increase performance I stored the calculated values and read them. This is faster than calculating all absolute values each time the program runs

# Storing the calculated data
  # meting1 <- LoadXsenseData(pp_info[1:3,1])
  # write.csv(meting1, file = "example1.csv", row.names = FALSE)

# Loading the calulated data
meting1 <- read.csv("example1.csv")

# Plot the calculated data
plot_a(meting1)

This look absolutely great! that the absolute acceleration approaches gravitational constant very well! It is a tiny bit higher than the expected 9.8 due to the noise. Since the absolute is taken, all noise that is not in opposite direction of the gravity increases the measured acceleration at rest. I’m overwhelmed by the precision of the IMU’s here!

Data clipping

The protocol was such that the barbell IMU only moved while the participant was executing a jerk, or while the barbell was being loaded. Thus, the barbell IMU data seems a wise place to identify the moments at which the jerk was executed.

Letst start by calculating the range of the baseline acceleration that we measured between t = 900 and t = 1100

# Find and print the minimum and maximum values from the range t = 900 - 1100
cat("minimum value of A_abs: ", min(meting1[meting1$time >= 900 & meting1$time <= 1100, ]))

minimum value of A_abs:  0.02195406

cat("maximum value of A_abs: ",max(meting1[meting1$time >= 900 & meting1$time <= 1100, ]))

maximum value of A_abs:  1100

Lets increase the interval slightly so that we can use it to automatically determine where movement occurs. We could exclude all cases where the barbell IMU measures an A_abs within the 8.5 - 11.2 inteval.However, when the barbell accelerates or decelerates the A_abs crosses the interval. Lets be conservative and assume that if it crosses the the interval, it does so for less than 360 elements (6 seconds! highly conservative but it works just fine!)

# New variable to work with
meting_filtered <- meting1

# Compute the run-length encoding
rle_sequence <- rle(meting_filtered$A_abs1 >= 8.5 & meting_filtered$A_abs1 <= 11.2)

# Identify the start and end indices of consecutive sequences where condition is TRUE
start_indices <- cumsum(rle_sequence$lengths) - rle_sequence$lengths + 1
end_indices <- cumsum(rle_sequence$lengths)

# Identify the consecutive sequences where the condition holds for more than 360 rows
condition_indices <- which(rle_sequence$values & rle_sequence$lengths >= 360)

# Iterate over the consecutive sequences and replace the values of A_abs2 with NA
for (i in condition_indices) {
  start_index <- start_indices[i]
  end_index <- end_indices[i]
  meting_filtered$A_abs1[start_index:end_index] <- NA
  meting_filtered$A_abs2[start_index:end_index] <- NA
  meting_filtered$A_abs3[start_index:end_index] <- NA
  meting_filtered$A_abs4[start_index:end_index] <- NA
  meting_filtered$A_abs5[start_index:end_index] <- NA
}

plot_a(meting_filtered)

rm(meting_filtered, rle_sequence, start_indices, end_indices, start_index, end_index, condition_indices, i)

The barbell IMU also registers acceleration while it is loaded. Since the subject was asked to sit still when resting, and the bar was only loaded while resting, it seems safe to assume that the IMU on the lower leg remained stationary (but with a more variable baseline) when no lift was exercised.

After visual inspection the interval of 8.5 - 11.2 seemed to suffice, again under the assumption of 360 elements

process_data1 <- function(meting_filtered) {
  # All columns where the conditions as described are true are removed. However, time is unchanged, since replancements work on indices. If an entire colum where to be removed, this would cause errors when plotting the data
  
  # Calculate the run-length encoding 
  rle_sequence_A_abs1 <- rle(meting_filtered$A_abs1 >= 9.2 & meting_filtered$A_abs1 <= 10.3)
  rle_sequence_A_abs2 <- rle(meting_filtered$A_abs2 >= 9.0 & meting_filtered$A_abs2 <= 11.5)
  rle_sequence_A_abs3 <- rle(meting_filtered$A_abs3 >= 9.0 & meting_filtered$A_abs3 <= 11.5)
  rle_sequence_A_abs4 <- rle(meting_filtered$A_abs4 >= 8.0 & meting_filtered$A_abs4 <= 11.5)
  rle_sequence_A_abs5 <- rle(meting_filtered$A_abs5 >= 7.2 & meting_filtered$A_abs5 <= 15.7)
  
  # Identify the start and end indices
  start_indices_A_abs1 <- cumsum(rle_sequence_A_abs1$lengths) - rle_sequence_A_abs1$lengths + 1
  end_indices_A_abs1 <- cumsum(rle_sequence_A_abs1$lengths)
  
  start_indices_A_abs2 <- cumsum(rle_sequence_A_abs2$lengths) - rle_sequence_A_abs2$lengths + 1
  end_indices_A_abs2 <- cumsum(rle_sequence_A_abs2$lengths)
  
  start_indices_A_abs3 <- cumsum(rle_sequence_A_abs3$lengths) - rle_sequence_A_abs3$lengths + 1
  end_indices_A_abs3 <- cumsum(rle_sequence_A_abs3$lengths)
  
  start_indices_A_abs4 <- cumsum(rle_sequence_A_abs4$lengths) - rle_sequence_A_abs4$lengths + 1
  end_indices_A_abs4 <- cumsum(rle_sequence_A_abs4$lengths)
  
  start_indices_A_abs5 <- cumsum(rle_sequence_A_abs5$lengths) - rle_sequence_A_abs5$lengths + 1
  end_indices_A_abs5 <- cumsum(rle_sequence_A_abs5$lengths)
  
  
  
  # Identify more than 10 consecutive indices
  condition_indices_A_abs1 <- which(rle_sequence_A_abs1$values & rle_sequence_A_abs1$lengths >= 300)
  condition_indices_A_abs2 <- which(rle_sequence_A_abs2$values & rle_sequence_A_abs2$lengths >= 360)
  condition_indices_A_abs3 <- which(rle_sequence_A_abs3$values & rle_sequence_A_abs3$lengths >= 360)
  condition_indices_A_abs4 <- which(rle_sequence_A_abs4$values & rle_sequence_A_abs4$lengths >= 360)
  condition_indices_A_abs5 <- which(rle_sequence_A_abs5$values & rle_sequence_A_abs5$lengths >= 900)
  
  # Merge the condition indices
  condition_indices <- unique(c(condition_indices_A_abs1, condition_indices_A_abs3, condition_indices_A_abs5))
  
  # Replace with NA if condition A_abs1 = true
  for (i in condition_indices_A_abs1) {
    start_index <- start_indices_A_abs1[i]
    end_index <- end_indices_A_abs1[i]
    meting_filtered[start_index:end_index, -1] <- NA  # -1 to exclude the 'time' column
  }
  
    # Replace with NA if condition A_abs2 = true
  for (i in condition_indices_A_abs2) {
    start_index <- start_indices_A_abs2[i]
    end_index <- end_indices_A_abs2[i]
    meting_filtered[start_index:end_index, -1] <- NA  # -1 to exclude the 'time' column
  }
  
  # Replace with NA if condition A_abs3 = true
  for (i in condition_indices_A_abs3) {
    start_index <- start_indices_A_abs3[i]
    end_index <- end_indices_A_abs3[i]
    meting_filtered[start_index:end_index, -1] <- NA  # -1 to exclude the 'time' column
  }
  
  # Replace with NA if condition A_abs4 = true
  for (i in condition_indices_A_abs4) {
    start_index <- start_indices_A_abs4[i]
    end_index <- end_indices_A_abs4[i]
    meting_filtered[start_index:end_index, -1] <- NA  # -1 to exclude the 'time' column
  }
  
  # Replace with NA if condition A_abs5 = true
  for (i in condition_indices_A_abs5) {
    start_index <- start_indices_A_abs5[i]
    end_index <- end_indices_A_abs5[i]
    meting_filtered[start_index:end_index, -1] <- NA  # -1 to exclude the 'time' column
  }
  return(meting_filtered)
}

# Call the function with your dataframe as argument
filtered1 <- process_data1(meting1)

plot_a(filtered1)

This distinguishes pretty well between lifts, except for the artifact at t=490, which we can remove manually

We now remain with some peaks that are due to the bar dropping, and this data should be excluded. As can bee seen, the barbells A_abs1 never really goes above ~20 when the hands are also moving (see for instance time 180 - 210). Obviously, the barbell should only be accelerating while the arms are also accelerating. Thus, the hig

filtered2 <- filtered1
filtered2 <- filtered2 %>%
  mutate(across(-time, ~ ifelse(time > 20 & time <= 30, NA, .))) %>%
  mutate(across(-time, ~ ifelse(time > 485 & time <= 530, NA, .))) %>%
  mutate(across(-time, ~ ifelse(time > 780, NA, .)))

plot_a(filtered2)

There remain two great peaks in the last two successful lifts. The subject dropped the bar and it hit the ground, where all the spead almost instanteniously dissapears, and thus A_abs1 peaks.Actually, the lift ends when the barbell is overhead, and thus the recording should stop here, regardless if the bar is dropped or delegated to front-rack. Therefore the second peak should somehow be identified and filtered out.

The arms lower when the bar is dropped or delegated to front-rack, but the peak acceleration of the arms is before the moment the bar hits the ground or the shoulders, and thus the onset of the final peak acceleration of the arms indicates the latest moment that the barbell was successfully overhead. Therefore, we can filter the data that comes after the onset of the final peak out.

After visual inspection of the data the arbitrary value of 60 was selected to determine movement of the arms. The final peak where A_abs5 reaches 60 is filtered out. peaks that are within 0.8 seconds of each other are regarded as one, and the 0.5 seconds before the peak are also discarded

# Additional package required
# 
# library(zoo)
# 
# # Define a function to find peaks
# find_peaks <- function(x) {
#   # Find peaks where A_abs5 is less than 60
#   peak_indices <- which(diff(sign(diff(x))) == -2) + 1
# 
#   # Remove peaks that are within 0.8 seconds of each other
#   peak_indices <- peak_indices[diff(df$time[peak_indices]) > 0.8]
# 
#   return(peak_indices)
# }
# 
# # Find peaks in A_abs5
# peaks <- find_peaks(df$A_abs5)
# 
# # Find the second peak
# second_peak <- peaks[2]
# 
# # Get the time of the second peak
# time_of_second_peak <- df$time[second_peak]
# 
# # Print the time of the second peak
# print(time_of_second_peak)